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MBE  Reconstruction  Using  p  =  400. 


I .  INTRODUCTION 


It  is  widely  appreciated  that  the  faithfulness  in  estimation  of  an 
unknown  object  depends  strongly  upon  what  is  known  a  priori  about  the 
object.  In  particular,  the  knowledge  that  the  unknown  object  must  be 
positive,  called  "positivity,"  has  been  well  exploited  in  restoration 
schemes. 1-10  Calculational  effort  with  reconstruction  techniques  is  also 


*A.C.  Schell,  "Enhancing  the  Angular  Resolution  of  Incoherent  Sources."  Radio 
Electronic,  Eng.,  Vol.  29,  pp.  21-26,  1965. 

P.A.  Jansson,  R.H.  Hunt,  and  E.K.  Plyler,  "Response  Function  for  Spectral 
Resolution  Enhancement,"  J.  Opt.  Soc.  Am.,  Vol.  58,  pp.  1665-1666,  1968. 

o 

Y.  Biraud,  "A  New  Approach  for  Increasing  the  Resolving  Power  of  Data 
Processing,"  Astron.  Astrophys.,  Vol.  1,  pp.  124-127,  1969. 

^B.R.  Friedon,  "Restoring  with  Maximum  Likelihood  and  Maximum  Entropy,"  J. 
Opt.  Soc.  Am.,  Vol.  62,  pp.  511-517,  1972. 

~*W.H.  Richardson,  "Bayesian  Based  Iterative  Method  of  Image  Restoration,"  J , 
Opt.  Soc.  Am.,  Vol.  62,  pp.  55-59,  1972. 

^L.B.  Lucy,  "An  Iterative  Technique  for  the  Rectification  of  Observed 
Distributions,"  Astron.  J.,  Vol.  79,  pp.  745-754,  1974. 

^S.J.  Wernecke  and  L.R.  D' Addario,  "Maximum  Entropy  Image  Reconstruction,” 
IEEE  Trans.  Computers,  Vol.  C-26,  pp.  351-364,  1977. 

Q 

S.F.  Gull  and  G.J.  Daniell,  "Image  Reconstruction  from  Incomplete  and  Noisy 
Data,"  Nature,  Vol.  272,  pp.  686-690,  1978. 

9 

B.R.  Frieden,  "Image  Restoration  Using  a  Norm  of  Maximum  Information,"  Proc. 
SPIE,  Vol.  207,  pp.  14-25,  1979. 

*®A.R.  Davies,  T.  Cochrane,  and  O.M.  Al-Faour,  "The  Numerical  Inversion  of 
Truncated  Autocorrelation  Functions,”  Optica  Acta,  Vol.  27,  pp.  107-118, 
1980. 
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considerably  reduced  when  non-physical  states  are  disallowed.  Positivity 

permits  one  to  produce  sharper  edge  gradients  where  the  edge  meets  the  known 
(or  zero)  background  level*”  whereby  the  estimated  gradient  profile  contains 
spatial  frequencies  that  may  appreciably  exceed  the  cutoff  frequency  in  the 
data.3-4 


However,  the  high  Intensities  near  the  top  of  the  edge  profile  are  not  so 
enhanced.  There,  the  data  are  far  from  the  "forbidden”  region  of  negative 
values;  enforcing  positivity  upon  already  positive  numbers  does  nothing  to 
them.  The  question  is,  then,  how  can  the  gradient  near  the  top  of  the  edge 
profile  be  enhanced  by  enforcing  some  other  prior  knowledge? 

Because  a  lower  bound  of  zero  works  at  low  intensity  values,  it  might  be 
expected  that  knowledge  of  a  finite  upper  bound  b  can  produce  the  desired 
effect  at  the  higher  intensity  values.  In  fact,  this  can  be  shown  to  be  true, 
merely  by  reapplying  the  argument  in  Reference  16  to  the  upper-bound 
situation.  The  enhancement  can  be  expected  to  be  most  effective  if  one  has  a 
least  upper  bound,  and  in  particular,  cases  where  the  object  attains  it  quite 
often  across  the  scene.  Ideally,  it  should  attain  the  bound  roughly  as  often 
as  it  attains  zero  or  background  values. 


Objects  of  this  kind  arise  diversely.  In  image  restoration,  the  object 
might  be  a  handwritten  or  printed  page,  where  the  white  paper  provides  the 
upper  bound  and  the  print  provides  the  lower  bound.  In  spectroscopy,  the 
object  can  be  an  absorption  spectrum,  which  must  lie  between  0  and  100% 
levels.*  In  photography,  advantage  be  taken  of  the  known  fog  and 
saturation  density  levels.*®  Or,  in  reconstruction  tomography,  the  object 


*  * R.  Gordon,  R.  Bender,  and  G.T.  Herman,  "Algebraic  Reconstruction  Techniques 
(ART)  for  Three  Dimensional  Electron  Microscopy  and  X-Ray  Photography,"  J. 
Theor.  Biol.,  Vol.  29,  pp.  471-481,  1970. 

*2G.  Minerbo,  "MENT,  A  Maximum  Entropy  Algorithm  for  Reconstructing  a  Source 
from  Projection  Data,"  Comp.  Graphics  Image  Processing,  Vol.  10,  pp.  48-68, 
1979.  ~~  . 

1 1 

A.  Lent,  "A  Convergent  Algorithm  for  Maximum  Entropy  Image  Restoration  with 
Medical  X-Ray  Applications,"  in  SPSE  Conference  Proceedings,  R.  Shaw,  ed., 
(SPSE,  Washington,  D.C.,  1976),  pp.  249-257. 

*4G.T.  Herman,  A.  Lent,  and  S.W.  Roland,  "ART:  Mathematics  and 
Applications,"  J.  Theor,  Biol.,  Vol.  42,  pp.  1-32,  1973. 

*3C.F.  Barton,  "Computerized  Axial  Tomography  for  Neutron  Radiography  of 
Nuclear  Fuel,”  Trans.  Amer.  Nucl.  Soc.,  Vol.  27,  pp.  212-213,  1977. 

^B.R,  Frieden,  "Estimation-A  New  Role  for  Maximum  Entropy,"  in  SPSE 

Conference  Proceedings,  R.  Shaw,  ed.  (SPSE,  Washington,  D.C.,  1976),  pp. 
261-265. 

*7P.A.  Jansson,  R.H.  Hunt,  and  E.K.  Plyler,  "Resolution  Enhancement  of 
Spectra,"  J.  Opt.  Soc.  Am.,  Vol.  60,  pp.  596-599,  1970. 

18  M 

B. R.  Frieden,  "Statistical  Estimates  of  Bounded  Optical  Scenes  by  the  Method 

of  Prior  Probabilities,"  IEEE  Trans.  Inform.  Theory,  Vol.  IT-19,  pp.  118- 
119,  1973.  “ 


might  be  machine  parts  or  rods  of  known  absorption  coefficient  immersed  in  a 
medium  whose  absorption  coefficient  is  also  known.  If  the  rods  are  denser 
than  the  medium,  they  provide  the  upper  bound  in  absorption,  the  medium 
providing  the  lower.  We  shall  simulate  this  case  in  particular  examples 
below. 


In  applications  of  this  type  one  may  have  to  choose  the  bounds  depending 
on  the  purpose  for  forming  the  tomographic  image.  If  the  purpose  is  to 
visualize  imperfections  in  the  rods,  it  is  important  that  the  imperfections  as 
well  as  the  rods  have  absorption  values  that  lie  between  the  prescribed  bounds 
(0,b).  Otherwise  the  Imperfections  will  not  be  augmented  by  enforcing  the 
bounds.  The  simplest  of  imperfections — cracks — satisfy  the  bounds  for  the 
rods  since  where  they  occur  a  level  b  has  merely  been  replaced  by  a  0.  This 
case  in  particular  will  be  considered  below  in  the  simulations.  If,  on  the 
other  hand,  the  imperfections  are  in  the  form  of  impurities  with  higher 
absorption  embedded  in  the  rods,  e.g.,  an  admixture  of  another  substance  that 
was  inadvertently  added  during  manufacture,  then  level  b  should  be  replaced  by 
a  proper  b'  exceeding  b. 

We  follow  the  approach  of  Reference  4  and  model  the  object  as  an  array  of 
photon  counts  n^,  n*l,  ...,  M=m,  where  mn  denotes  the  number  of  photons 
absorbed  (in  the  tomographic  case  of  interest)  at  position  xn  in  the  object. 
Let  -  o  represent  the  energy  increment  represented  by  each  count.  Then  an 
object  o  relates  to  its  counts  through 

o  =  mAo.  (la) 

We  seek  the  most  probable  object  o  consistent  with  two  known  bounds 

0  <  o^<  b,  n  =  l,  • • • ,M.  (lb) 

II.  IMAGE  MODELING 


A.  Object  Model 

In  this  section  we  discuss  a  model  of  the  object.  Let  the  photons  in 
question  be  x  rays,  as  for  example  in  tomography.  These  have  low  quantum 
degeneracy^  and  hence  behave  statistically  like  discrete  or  Boltzmann 
particles.  This  is  one  aspect  of  the  model.  Another  aspect  is  the  mean 
object  <o> .  By  the  law  of  large  numbers,^®  as  the  number  of  photons 
approaches  infinity  o  will  approach  <o>.  Hence  the  numbers  <o>  act  as 
"biases"  for  o.  Since  biasing  is  a  profound  effect,  it  must  be  closely 


19 

R.  Kikuchi  and  B.H.  Soffer,  "Maximum  Entropy  Image  Restoration.  I.  The 
Entropy  Expression,"  J.  Opt.  Soc.  Am..  Vol.  67,  pp.  1656-1665,  1977. 

20 

B.R.  Frieden,  Probability,  Statistical  Optics  and  Data  Testing, 
Springer-Verlag,  New  York,  1983. 
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considered.  If  a  form  is  assumed  for  <o> ,  then  o  can  only  fluctuate 
probabilistically  about  it.  How  can  we  model  <o>? 


An  important  aspect  of  Che  estimate  o  will  be  its  reliability.  Again, 
using  the  example  of  rods,  if  o  shows  a  crack  or  pit,  can  it  be  trusted?  The 
detail  called  "a  crack”  consists  of  a  strong  local  departure  from  grayness 
that  is,  it  has  an  abrupt  decrease  from  b  to  0  and  corresponding  large 
gradient.  Hence,  if  we  knew  that  o  only  shows  small  departures  from  grayness, 
we  could  trust  the  crack  detail.  On  the  other  hand,  we  found  that  numbers  <o> 
act  as  biases  for  o.  It  follows  then  that  if  numbers  <o>  were  flat  or 
constant 


<on>  -  y  =•  constant. 


and  if  we  knew  this  to  be  true,  then  reliability  could  be  so  built  into  the 
estimated  o. 


The  requirement  (2)  can  be  nearly  true  if  the  object  consists  of  half 
background  (with  but  a  few  rods),  but  it  will  not  be  rigorously  true  unless 
the  object  is  only  background  (a  most  uninteresting  case).  However,  the 
approximation  has  been  used  successfully  to  process  astronomical  pictures. 


Similarly,  if  the  object  field  is  packed  with  rods,  requirement  (2)  will 
be  far  from  satisfied.  In  this  case  art  assumption  of  (2)  will  produce  errors 
in  the  estimate;  however,  these  will  tend  to  be  errors  of  "omission"  only. 
That  is,  certain  details  will  be  missed,  but  no  artifacts  will  tend  to  be 
created.  Again,  the  reconstruction  of  a  crack  will  be  viewed  as  probably 
truthful.  However,  some  cracks  might  be  missed  as  the  price  paid.  The  false 
alarm  rate  will  be  low,  but  perhaps  so  will  be  the  detection  rate.  Empirical 
tests  may  establish  the  rates  in  particular  applications.  (Note:  empirically 
the  theoretical  possibility  of  a  low  detection  rate  was  not  borne  out.  See 
the  Applications  section  below.)  We  shall  adapt  condition  (2)  as  the  second 
aspect  of  the  object  model. 


B.  A  Priori  Probability  of  an  Object 


Given  the  particle-like  behavior  of  the  photons,  the  knowledge  of  bounds 
(lb)  and  the  assumption  of  a  gray  mean  object  (2),  we  are  ready  to  form  Pl(o), 
the  a  priori  probability  of  a  photon-count  object  o.  This  will  be  the 
probability  that  m  indistinguishable  particles  are  distributed  in  any  order 
within  M  cells,  where  each  cell  can  hold  anywhere  from  0  to  b/Ao  particles. 
This  obeys  the  binomial  statistic 


„  ft  (b/Ao)!  m  b/ Ao-m_ 

P  (o)  =  II  — mm - VT  P  n  q  n 

I  n=l  in  !  b/Ao  -  m  !  rn  n  , 


m  =  ■*—  , 

n  Ao 


q  =  1  -  p  . 
'n  '  n 


^B.R.  Frteden  and  D.C.  WelLs,  "Restoring  with  Maximum  Entropy.  III.  Poisson 
Sources  and  Backgrounds,”  .1 .  Opt .  Soc .  Am.  ,  Vol.  68,  pp.  93-103,  1978. 


10 


In  effect,  the  x  rays  are  being  treated  like  electrons,  since  this  is  also  the 
likelihood  expression  for  Fermi-Dirac  particles.  [Note:  A  more  rigorous 
derivation  of  (3)  from  the  standpoint  of  photon  (not  particle)  statistics  is 
given  in  Appendix  A.]  In  particular  for  the  mean  object  (2),  this  becomes 


P.eo  -  a/2>Mb/4°  »  , 
l  n=  1 


(b/Ao) ! 

m  ! (b/ Ao  -  m  ) ! 
n  n 


const,  x  (o  /Ao)!(b/Ao  -  o  /Ao)!  * 
n  n 


the  sought  expression. 


C.  Image  Model 


The  rest  of  the  theory  in  this  paper  follows  that  of  a  previous 
paper. ^  The  theory  makes  the  following  assumptions  about  the  image  data 
formed  from  the  unknown  object  o: 


(1)  The  image  i,,  ...,  i^=i  suffers  from  noise  , 
known  background  profile  Bj  ,  . . . ,  15  such  that 

i  =  ?,  o  s(x  -  x  )  +  B  +  n  . 
m  n=l  n  m  n  ra  m 


Oyp=n  and  rides  atop  a 


The  quantity  s  is  the  point  spread  function  of  the  imagery.  In  computer 
tomography,  it  is  the  rayed  function  resembling  the  British  Union  Jack,  with 
one  ray  for  each  projection. 

Thus,  there  are  now  two  sets  of  unknowns,  o  and  n.  Background  B  is 
assumed  known  by  the  use  of  some  prefiltering  operation  upon  the  image,  such 
as  median  windowing  it.^-24  (2)  The  noise  n  is  Poisson.  This  is  indeed  the 

case  when  using  modern  imaging  arrays  in  astronomy^  and  in  computer 
tomography. ^ 


2  y  M 

B.R.  Frieden,  "New  Restoring  Algorithm  for  the  Preferential  Enhancement  of 
Edge  Gradients,"  J,  Opt.  Soc.  Am.,  Vol.  66,  pp.  280-283,  1976. 

21 

N.C.  Gallagher  and  G.L.  Wise,  "Passband  and  Stopband  Properties  of 
Median  Filters,”  Proceedings  of  the  1980  Conference  on  Information  Sciences 
and  Systems,  Princeton  University,  New  Jersey,  1980,  pp.  303-307. 

^B.R.  Frieden,  "Some  Statistical  Properties  of  the  Median  Window,"  Proc , 
SPIE,  Vol.  373,  pp.  2 19-224 ,  1981. 

2 

H.H.  Barrett  and  W.  Swindell,  Radiological  Imaging,  Academic  Press, 

New  York,  1981. 
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III.  THE  ALGORITHM 


A.  Net  Likelihood  Function 


Taking  a  conventional  probabilistic  approach,  we  seek  unknowns  o  and  n 
that  are  jointly  maximum  probable, 


P(o,n)  =  maximum, 


(6) 


subject  to  the  data.  From  elementary  considerations 


P(o,n)  =  P1(o)P2(n|o) . 


C7) 


We  already  know  P^(o);  see  Equation  4. 


The  conditional  probability  P2(i|o)  defines  the  fluctuations  in  i  given 
one  object  o.  The  noise  was  assumed  to  be  Poisson  on  the  image  with  the  image 
given  as  photon  counts.  If  each  count  consists  of  an  energy  increment  i, 
then  the  image  intensity  i  corresponds  to  i/ Ai  counts.  Then,  assuming 
independent  image  values,  we  have 


i  /Ai  -a 
„  a  m  e  m 

P2(1I0)  =  tJl  (i  /Ai)! 

m 


=  Pinlo) 
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where  a  is  the  noiseless  signal  image  count.  By  Equation  (5), 
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The  final  identity  in  (8)  follows  because  with  o  fixed,  by  Equation  (5) 
corresponding  values  of  im  and  n  have  the  same  histogram.  The  net  likelihood 
function  is  then,  by  identity  (7),  the  product  of  Equations  (4)  and  (8). 

B.  Restoring  Principle 

The  principle  of  restoration  is  to  maximize  P(o,n)  through  choice  of  o 
and  n  subject  to  the  image  data  idata  obeying  (5).  We  shall  also  assume  that 
the  total  energy  E  in  the  object  is  known,  e.g.,  by  conservation  of  energy 
from  the  image  data.  It  is  mathematically  convenient  to  maximize  In  P(o,n), 
instead  of  P(o,n)  which  gives  the  same  solution.  Also,  we  add  the  data  and 
energy  constraints  to  the  objective  function  via  Lagrange 

multipliers  X  and  y  .  Accordingly,  by  Equations  6  and  7  we  have  to  maximize 
the  function 


In  Pj(o)  +  In  P(n  o) 


1  X  (i  -  idata)  -  p(5;  o  -  E) 
m=  1  m  m  m  n 


ao) 


through  choice  of  n  ,o  and  the  Lagrange  multipliers.  Then  by  Equations  (4) 
and  (8)  the  objective  function  is 
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(11) 


£  (o  /Ao)  ln(o  /ho)  -  E  (b/Ao  -  o  / ho)  ln(b/Ao  -  o  /ho) 
n  n  n  n  n  n 


+  In  P(n|o)  -  E  X  (i  -  i  ata)  -  p(I  o  -E)=*  maximum  . 

1  m  m  m  m  n 

The  first  sum  is  the  Shannon  entropy  of  the  object,  while  the  second  is  the 
Shannon  entropy  of  the  unfilled  photon  sites  within  each  pixel.  In  other 
words,  we  have  filled  entropy  plus  unfilled  entropy.  This  appears  to  be  the 
natural  way  to  accommodate  an  upper  bound  into  a  maximum  entropy  approach. 
(See  also  Reference  18.) 


C.  Net  Restoring  Algorithm 

The  solution  to  (11)  is  f ound^  ^  by  substituting  in  the  expression  (11) 
Equation  (8)  for  P(n|o)  and  Equation  (5)  for  im;  and  by  the  usual  rules  of 
calculus,  equating  to  0  in  turn  the  partial  derivatives  of  the  objective 
function  3/ So^  (n  fixed)  and  3/ 3n  (o  fixed).  Also,  two  approximations  are 
made:  a  large  enough  number  of  photons  are  assumed  present  so  that  the 
Poisson  law  (8)  may  be  approximately  replaced  by  a  normal  law  whose  variance 
equals  the  mean;  and  the  signal  image  a  is  assumed  to  be  smooth  and  slowly 
varying.  With  these  approximations  the  solution  is  an  o  obeying 


jdata  _  [  jM  o  s(x  -  x  )  +  B  ]  e  *  ^n/P,  m  =  1 ,  . . .  ,M  (12a) 

m  n=l  n  m  n  m 

E  =  Jjj  on,  (12b) 

o  =  - - -  ,  n  =  1,  ...,  M  (12c) 

1  +  exp  [T+EAs(x  -x)] 
m  m  m  n 


r  =  uAo,  A  =  X  Ao,  p  =  Ao/Ai.  (12d) 

’mm’ 

Thus ,  o  represented  through  (12c)  in  terms  of  M  +  1  free  parameters  T,  A,  must 

obey  the  M  +  1  data  Equations  12a, b.  This  is  our  maximum  bounded  entropy 

( M BE)  restoring  algorithm.  We  can  see  from  the  form  of  Equation  12c  that  the 
estimated  on  cannot  have  values  outside  the  interval  (0,b).  The  prescribed 
quantity  p,  defined  at  (12d),  is  called  the  "sharpness  parameter."  A  higher 
value  of  p  increases  the  resolution  of  the  output  o.  Typically  a  value  p  **  50 

causes  a  modest  increase  in  resolution,  while  p  =  200  causes  a  high 

increase.  This  behavior  is  consistent  with  its  definition  in  ( 1 2d ) :  If  one 
inputs  a  high  p  ,  he  is  assuming  that  the  object  intrinsically  consists  of 
jumps  Ao  in  intensity  much  exceeding  those  in  the  given  image.  This  causes  a 
"jumpier"  and  hence  higher-resolved  output.  Other  properties  of  p  are 
discussed  in  Reference  21. 


As  mentioned  before,  B  is  the  estimated  background  intensity  function. 
Background  is  defined  as  a  slowly  varying  component  of  the  image  data  which 
intrinsically  lacks  resolution  and  hence  is  incapable  of  being  restored 
further.  The  function  B  may  be  estimated  by  purposely  blurring  the  input 
image,  either  by  convolution  with  (say)  a  pillbox  function^  or  by  the  use  of 
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a  median  window.  The  latter  approach  is  recommended  when  the  object  details 
of  interest  have  a  known  largest  support.  In  this  case,  the  use  of  a 
circular,  filled  •  Jndow  of  diameter  equal  to  twice  (or  more)  this  support 
value  should  be  used.  The  output  of  this  operation  is  "blind"  to  these  object 
details,  and  hence  only  "sees"  the  background.  Although  this  approach  is  more 
time-consuming  than  the  convolution  approach  mentioned  above,  it  gives  a  more 
accurate  estimate  of  the  background. 

IV.  APPLICATIONS  TO  COMPUTER  TOMOGRAPHY 

The  preceding  algorithm  has  been  developed  for  any  imaging  situation 
where  (a)  the  photons  behave  like  particles,  (b)  the  object  is  bounded  by 
intensity  levels  0  and  b  where  b  is  a  least  upper  bound  to  intensity,  (c)  the 
object  consists  of  a  slowly  varying  background  function  plus  a  foreground 
function  whose  details  it  is  desired  to  restore,  (d)  the  image  is  formed 
convolutionally  from  the  object  via  a  known  point  spread  function,  and  (e)  the 
image  suffers  from  Poisson  noise. 

There  are  many  cases  where  these  conditions  are  satisfied  in  particular 
the  case  of  computer  tomographic  imaging.  Suppose  that  the  "back-projected" 
linage^  is  the  given  data  i^ata.  This  is  known  to  connect  with  the 
absorptance  object  o  via  convolution  with  a  Union  Jack  point  spread  function, 
for  example  such  as  in  Figure  la.  Ea«_..  arm  of  the  pattern  is  formed  by  a 
projection  in  that  direction.  Hence  in  this  example  we  are  working  with  four 
projections. 

The  objects  of  interest  are  rods  immersed  in  a  medium.  The  rod  cross 
sections  comprise  our  foreground  details.  The  absorptances  of  the  rods  are 
known  to  be  at  level  b,  i.e.,  all  are  at  the  upper  bound;  clearly,  this  is  an 
ideal  application  of  the  approach.  The  absorptance  of  the  medium  is  also 
known,  but  for  consistency  in  the  approach  we  estimate  it  (below).  A  typical 
object  of  this  type  is  shown  in  Figure  lb,  where  the  object  rods  are  shown 
within  a  large  cylinder.  Everything  beyond  the  cylinder  is  at  0  level. 
Intensity  levels  have  been  logarithmically  stretched  so  as  to  enable  the 
background  to  be  seen:  it  is  at  5%  of  the  foreground.  Hence  the  object  has 
high  contrast. 

Notice  that  the  foreground  rod  cross  sections  consist  of  three  shapes 
approximating  circles,  the  largest  on  top,  smallest  to  the  lower  left.  Each 
pixel  is  one  detector-width  wide,  i.e.,  contains  only  one  ray  from  a  given 
projection.  Hence,  the  back-projected  image  of  this  object  will  suffer  severe 
spillover  of  energy  from  the  broadest  rod  into  the  other  two,  and  vice 
versa.  In  other  words,  there  will  be  severe  blur  present.  This  was  to 
provide  an  "acid-test"  of  the  approach. 

The  spread  function  1(a)  was  convolved  in  the  computer  with  1(b)  to 
produce  the  back-projected  image.  This  was  made  a  Poisson  noise  process  with 
a  signal-to-noise  ratio  (S/N)  of  10:1  at  the  brightest  pixel  in  the  image, 
Figure  lc.  The  image  is  quite  noisy,  since  S/N  falls  off  as  the  square  root 
of  intensity  and  hence  is  much  less  than  10:1  at  most  points  in  the  scene. 

The  principal  blur  is  visually  along  the  four  projection  directions,  as  would 
be  expected.  This  image  also  suffers  from  numerous  artifact  "sources"  due  to 
the  chance  crossing  of  rays  from  different  projections. 


14 


To  speed  up  the  MBE  algorithm  (12a-d)  we  programmed  it  assuming  a 
Gaussian  point  spread  function  s  to  be  present.  A  Gaussian  spread  function  is 
separable  in  the  x  and  y  directions.  This  permits  the  algorithm  to  be  applied 
one  dimensionally,  first  to  the  rows,  then  to  the  columns,  of  the  image.21 
But  the  actual  point  spread  function  is  not  Gaussian,  as  is  seen  in  Figure 
la.  Hence,  we  had  to  filter  the  Poisson  image  into  a  Gaussian  form, 
essentially  dividing  out  the  Union  Jack  and  multiplying  by  the  Gaussian.  This 
output  image  is  called  the  "prefiltered  image."  Using  a  standard  deviation  => 

2  pixels  in  the  Gaussian  spread  function,  we  obtained  Figure  Id.  Notice  that 
the  rayed  appearance  of  the  Poisson  image  1(c)  has  been  somewhat  reduced. 

This  is  due  to  the  Gaussian  falloff  of  the  new  point  spread  function. 

However,  some  artifacts  linger  on.  This  image  comprises  idata>  the  input  to 
the  MBE  algorithm. 

Since  the  image  has  been  filtered,  the  background  region  has  been  changed 
from  its  value  in  the  object.  Hence,  it  has  to  be  estimated.  To  do  this,  we 
took  advantage  of  knowing  that  the  rods  are  rather  packed  near  the  center,  so 
that  beyond  a  radius  of  about  16  pixels  only  background  occurs  (out  to  the 
cylinder  walls).  Hence,  an  average  was  taken  over  this  region  in  1(d)  to 
infer  the  new  background  level.  Background  function  B  was  then  made  to  be 
this  constant  value  out  to  the  cylinder  walls,  and  thereafter  the  image  1(d) 
itself;  the  latter  because,  beyond  the  v’-ita  there  is  known  to  be  no 
foreground  object.  This  background  image  is  shown  in  Figure  le.  Numerous 
artifacts  in  the  form  of  "blobs"  can  be  seen  outside  the  walls.  These  are  due 
to  imperfect  prefiltering  of  the  Union  Jack  into  the  Gaussian.  However,  they 
are  in  actuality  much  weaker  than  seen:  we  gray-scale  stretched  so  as  to 
render  them  visible.  Also,  since  these  lie  outside  the  region  of  interest 
(the  foreground  object  details),  they  do  not  much  interfere  with  the  MBE 
outputs. 

Knowing  i^3*-3  and  B,  we  can  now  use  the  MBE  routine.  The  use  of  a 
sharpness  parameter  200  resulted  in  the  restoration  shown  in  Figure  If.  This 
is  a  pretty  fair  reconstruction  of  the  original  object  1(b).  Its  good  aspects 
are  (a)  complete  resolution  of  the  three  rods;  compare  with  the  resolution 
present  in  1(c)  or  1(d);  (b)  an  almost  absence  of  artifacts  (one  is  visible  on 
the  lower  left);  (c)  very  strong  edge-gradients  at  the  rod  boundaries;  (d) 
true  absorptance  values  within  the  restored  rods  (white  corresponds  to  level 
b);  and  (e)  faithful  reconstruciton  of  the  top  rod's  shape,  most  probably 
because  it  has  the  most  energy  of  the  three  and  hence  suffered  from  noise 
propagation  the  least.  The  bad  aspects  are  (a)  faulty  shape  in  the 
reconstructed,  lower-right  rod;  and  (b)  strongly  underestimating  the  intensity 
in  the  lower-left  rod — it  is  almost  not  visible.  But,  considering  that  these 
results  followed  from  the  use  of  only  four  projection  directions,  we 
considered  them  encouraging. 

We  proceeded  to  try,  in  the  same  way,  the  case  of  20  projections. 
Corresponding  results  are  shown  in  Figure  2(a-e).  In  Figure  2e  is  shown  the 
MBE  output  for  p  =  200.  This  is  superior  to  the  corresponding  four- 
projections  output  in  Figure  If,  as  was  expected.  In  particular,  resolution 
is  very  high,  shapes  are  more  faithfully  reconstructed,  artifacts  are  still 
low,  and  even  the  weak,  lower-left  source  is  (now)  strongly  restored. 

We  also  compared  these  results  with  results  by  the  maximum  entropy  (ME) 
algorithm  for  the  same  data.  Notice  that  In  the  objective  function  (11)  if  b 


is  made  very  large  the  second  sum  is  effectively  independent  of  choice  o,  i.e. 
it  is  a  constant.  Therefore,  in  this  case  the  algorithm  simply  maximizes 
entropy,  the  first  sum  and  the  MBE  algorithm  becomes  the  ME  algorithm.  Hence, 
we  made  b  twice  the  known  upper  bound  (any  further  increase  did  not 
significantly  change  the  output).  The  resulting  ME  restoration  is  shown  in 
Figure  2f.  This  is  a  softer,  lower-resolution  estimate  than  before.  The 
lower-left  object  is  now  barely  detected.  We  conclude  that  the  MBE  algorithm 
with  a  good  estimate  of  the  upper  bound  b  has  a  strong  advantage  over  the 
maximum  entropy  (ME)  algorithm. 

In  order  to  observe  the  effect  of  increased  p  upon  the  reconstruction,  we 
now  increased  p  to  a  value  300,  with  output  in  Figure  2g,  and  to  a  value  400, 
with  output  in  Figure  2h  (level  b  was  once  more  at  its  true  value).  There  is 
some  increase  in  resolution,  but  it  is  apparent  that  this  is  saturating.  The 
envisioned  application  to  rod  cross  sections  is  to  detect  defects, 
particularly  cracks,  in  the  rods.  As  an  example  we  included  a  one  pixel  wide 
crack  in  a  rod.  So  narrow  a  crack  furnishes  an  acid  test  for  the  algorithm. 

We  again  used  20  projections.  Results  are  shown  in  Figure  3.  The  cracked  rod 
is  the  top  one,  as  shown  in  3b.  The  overall  object  is  otherwise  the  same  as 
in  Figure  2b.  The  Poisson  image  in  3(c)  does  show  a  gray,  nebulous  shape  in 
the  vicinity  of  the  crack.  However,  this  cannot  be  definitively  used  as  an 
indicator  for  a  crack,  since  image  2(c)  also  has  a  gray,  nebulous  shape  there 
whereas  its  object  did  not  have  a  cr.  r .  evidently,  some  of  the  gray  shape  is 
due  to  the  particular  noise  values  chosen,  which  are  about  the  same  in  both 
cases. 

We  prefiltered  the  Poisson  image  into  its  Gaussian  form  in  Figure  3d,  now 
using  a  o  of  1.5  (smaller  than  in  Figure  2d).  We  then  restored  this  by  MBE  in 
3(e),  using  a  p  of  200.  This  again  sharply  restores  the  rod  cross  sections, 
but  now  with  a  notch  in  the  top  rod.  This  notch  exactly  corresponds  in 
position  to  the  crack.  Comparison  with  the  corresponding  rod  reconstruction 
in  Figure  2e,  where  the  object  did  not  have  a  crack,  shows  a  decisive 
difference.  It  is  apparent  that  the  crack  has  definitely  been  reconstructed 
in  3(e) . 

Figure  3f  shows  the  result  of  now  using  MBE  with  p  =  400.  Higher 
resolution  is  attained,  with  the  crack  slightly  better  restored.  This  should 
be  compared  with  reconstruction  Figure  2h,  where  the  crack  did  not  exist  in 
the  object. 

A  better  comparison  with  results  in  Figure  2  is  obtained  if  the  same  a  as 
in  Figure  2  were  used  in  the  prefiltering  step  in  both  cases.  Accordingly,  we 
filtered  image  3(c)  into  a  Gaussian  form  using,  a  =  2.  The  result  is  Figure 
3g.  This  image  was  fed  into  MBE  using  p  =  400,  with  the  result  3(h). 
Comparison  with  the  corresponding  uncracked  restoration  Figure  2h  shows  a 
definite  restoration  of  the  crack,  although  not  quite  as  vividly  as  in  3(e)  or 
3(f).  Evidently,  the  use  of  a  smaller  o  helped  in  this  example. 

We  may  summarize  these  results  by  stating  that  MBE  can  restore  one  pixel 
wide  cracks  in  the  rods  with  high  reliability.  This  result  is  obtained  for 
the  high-contrast  objects  tested  here,  assuming  20  or  more  projections,  and 
with  S/N  the  order  of  10:1  (or  better).  Thus,  the  biasing  of  the  outputs 
toward  a  flat,  gray  scene  did  not  cause  the  algorithm  to  miss  crack  details 
under  these  conditions.  (See  Object  Model  section.)  We  have  not  tested  the 
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algorithm  on  low-contrast  objects  or  with  lower  S/N  in  the  image.  In  the 
tested  case,  the  crack  was  only  partially  restored,  producing  a  distorted 
image  of  the  rod.  In  practical  applications  it  might  not  be  recognized  as  a 
crack,  since  also  the  smaller  undamaged  rods  were  distorted.  However,  one  can 
expect  a  good  restoration  of  cracks  wider  than  one  pixel. 

V.  CONCLUSIONS 

The  knowledge  of  a  least  signal  upper  bound  is  very  effective  in 
enhancing  the  resolution  of  image  reconstructions.  The  proviso  is  that  this 
bound  be  met  over  a  substantial  part  of  the  object  field.  The  number  of 
projections  that  are  needed  for  high  quality  outputs  by  the  suggested  method 
can  be  quite  small,  ranging  from  4  to  20,  depending  on  accuracy 
requirements.  In  a  sample  problem  rod  cross  sections  could  be  accurately 
reconstructed  with  very  high  edge  gradients,  and  one-pixel-wide  cracks  can  be 
restored. 

The  time  requirements  for  the  64  x  64  pixel  cases  shown  were  about  8  s  of 
CPU  time  on  a  Cyber  135  mainframe  computer.  The  time  requirement  is 
proportional  to  the  area  of  the  image  to  be  processed  expressed  in  pixels. 
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APPENDIX  A:  USE  OF  BOSE-EINSTEIN  STATISTICS 


A  more  fundamental  derivation  of  the  a  priori  object  law  (3)  follows  by 
separately  considering  photon  absorptions  and  transmissions  through  the 
object .* 

Let  En  known  photons  be  incident  upon  pixel  n.  This  quantity  is 
ordinarily  known  in  computer  tomography,  since  the  incident  energy  upon  the 
object  is  accurately  monitored,  and  pixels  toward  the  emergent  side  of  the 
object  are  not  strongly  blocked  by  pixels  on  the  incident  side  (most  of  the 
light  passes  unabsorbed  through  the  object). 

Next,  consider  the  identity 

En  =  m  +  (b/o  -  m)  +  (Ejj  -  b/ Ao)  .  (Al) 

This  describes  the  fate  of  the  E  photons  as  independently  m  absorptions,  (b/o 
-  m)  transmissions  and  (En  -  b/o)  transmissions.  Let  these  three  states  have 
z,  z' ,  and  z”  degrees  of  freedom,  respectively.  Then  the  probability  law  for 
the  En  photons  is 

P(En)  =  Pz(m)Pz’(b/Ao  -  m)Pz"(En  -  b/ Ao)  ,  (A2) 

where  Pz(k)  is  proportional  to  the  Bose-Einstein  statistics,^ 

Pz(k)  =  ((k  +  z  -  l)!)/(k!(z  -  1)!)  pk.  (A3) 

By  inspection,  in  the  classical  particle  limit  z  large,  this  law  goes  over 
into  a  statistic 

Pz(k)  =  pk/k!  .  (A4) 

The  probability  law  Pz"  in  (A2)  is  merely  a  multiplicative  constant 
independent  of  m  and  here  may  be  ignored.  Hence,  in  the  classical  particle 
limit  z,  z'  large  appropriate  for  x  rays,  Equation  (A2)  becomes  proportional 
to  a  binomial  statistic, 

P(En)  a  pnm/m!  [qnb/o-ra/(b/  Ao  -  m) ] !  .  (A3) 

As  the  pixels  n  act  independently,  the  net  probability  P(E^,...,Em)  for  all 
photons  goes  over  into  a  product  of  factors  (A5).  This  is  of  the  form 
Equation  3,  as  was  to  be  proved. 


*We  thank  B.H.  Soffer  of  Hughes  Research  Laboratories  for  the  basic  idea 
behind  this  proof. 
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Cambridge,  MA  02139 

1  University  of  North  Carolina 

Bio-Medical  Engineering 
Department 

ATTN:  F. A.  DiBianca 
Chapel  Hill,  NC  27514 

1  Hospital  of  the  University  of 

Pennsylvania 
Department  of  Radiology 
ATTN:  G.T.  Herman 
3400  Spruce  Street 
Philadelphia,  PA  19104 

1  University  of  Pittsburgh 

ATTN:  D.  Sashin 
Pittsburgh,  PA  15213 

1  Stanford  University 

Department  of  Physics 
ATTN:  A.  Macovski 
Palo  Alto,  CA  94305 

1  University  of  Utah 

Department  of  Mathematics 

ATTN:  F.  Stenger 

Salt  Lake  City,  UT  84112 

1  Dr.  G.  Morgan 

1912  Sheffield  Court 
Severn,  MD  21144 

2  Mathematics  Research  Center 
ATTN:  B.  Noble 

J.  Nohel 

610  Walnut  Street 
Madison,  WI  53706 

1  Advanced  Research  and 

Applications  Corporation 
ATTN:  A.G.  Lason 
8150  Leesburg  Pike 
Vienna,  VA  22180 
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General  Electric  Armanent 
5  Electrical  Systems 
ATTN:  M.  J.  Bui man 
Lakeside  Avenue 
Burlington,  VT  05401 

Bell  Laboratories 
ATTN:  M.  Sondhi 
F.  Shepp 

Murray  Hill,  NJ  07971 
BDM  Corp. 

ATTN:  T.P.  Goddard 

2600  Cearden  Road 
North  Bldg 

Monterey,  CA  93940 

General  Electric  Computer 
Research 

Medical  Diagnostics  Systems 
Program 

ATTN:  T.  Kincaid 
Schenectady,  NY  12345 

Shell  Development  Oil  Co. 
ATTN:  Dr.  H.  Vinegar 
P.O.  Box  481 
Houston,  TX  77001 

Princeton  University 
Forrestal  Campus  Library 
ATTN:  K.  Brezinsky 

I.  Glassman 
P.O.  Box  710 
Princeton,  NJ  08540 

Princeton  University 
MAE  Dept. 

ATTN:  F.A.  Williams 
Princeton,  NJ  08544 

Purdue  University 
School  of  Aeronautics 
and  Astronautics 
ATTN:  R.  Click 

J. R.  Osborn 
Grissom  Hall 

West  Lafayette,  IN  47906 
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Purdue  University 
School  of  Mechanical 
Engineering 

ATTN:  N.M.  Laurendeau 
S.N.B.  Murthy 
TSPC  Chaffee  Hall 
West  Lafayette,  IN  47906 

Rensselaer  Polytechnic  Inst. 
Dept,  of  Chemical  Engineering 
ATTN:  A.  Fontijn 
Troy,  NY  12181 

Southwest  Research  Institute 
ATTN:  R.E.  White 

A. B.  Wenzel 
8500  Culebra  Road 
San  Antonio,  TX  78228 

Stanford  University 
Dept,  of  Mechanical 
Engineering 
ATTN:  R.  Hanson 

Stanford,  CA  94305 

University  of  Texas 
Dept,  of  Chemistry 
ATTN:  W.  Gardiner 
Austin,  TX  78712 

University  of  Utah 

Dept,  of  Chemical  Engineering 

ATTN:  G.  Flandro 

Salt  Lake  City,  UT  84112 

Virginia  Polytechnic 
Institute  and 
State  University 
ATTN:  J.A.  Schetz 
R.T.  Smith 

Blacksburg,  VA  24061 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the 
reports  it  publishes.  Your  comments/answers  to  the  items/questions  below  will 
aid  us  in  our  efforts. 

1 .  BRL  Report  Number _ Date  of  Report 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or 

other  area  of  interest  for  which  the  report  will  be  used.) _ 


4.  How  specifically,  is  the  report  being  used?  (Information  source,  design 
data,  procedure,  source  of  ideas,  etc.) _ 


S.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far 
as  man-hours  or  dollars  saved,  operating  costs  avoided  or  efficiencies  achieved 
etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future 
reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 


Name 


CURRENT 

ADDRESS 


Organi zat ion 


Address 


City,  State,  Zip 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the 
New  or  Correct  Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 
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('hD  Organization 

ADDRESS 

Address 


City,  State,  Zip 


(Remove  this  sheet  along  the  perforation,  fold  as  indicated,  staple  or  tape 
closed ,  and  mail . ) 


Director 

US  Army  Ballistic  Research  Laboratory 
ATTN:  AMXBR-OD-ST 

Aberdeen  Proving  Ground,  MD  21005-5066 


OFFICIAL  BUSINESS 

PENALTY  FOR  PRIVATE  USE.  *300 


BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  NO  12062  WASHINGTON, DC 
POSTAGE  WILL  BE  PAID  BY  DEPARTMENT  OF  THE  ARMY 


Director 

US  Army  Ballistic  Research  Laboratory 
ATTN:  AMXBR-OD-ST 

Aberdeen  Proving  Ground,  MD  21005-9989 
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